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We describe general characteristics of the Hodgkin-Huxley neuron's response to a periodic train 
of short current pulses with Gaussian noise. The deterministic neuron is bistable for antiresonant 
frequencies. When the stimuli arrive at the resonant frequency the hring rate is a continuous function 
of the current amplitude lo and scales as (/o — IthY^'^, characteristic of a saddle-node bifurcation 
at the threshold Ith- Intervals of continuous irregular response alternate with integer mode- locked 
regions with bistable excitation edge. There is an even-all multimodal transition between the 2:1 
and 3:1 states in the vicinity of the main resonance, which is analogous to the odd-all transition 
discovered earlier in the high-frequency regime. For Iq < Ith and small noise the firing rate has a 
maximum at the resonant frequency. For larger noise and subthreshold stimulation the maximum 
firing rate initially shifts toward lower frequencies, then returns to higher frequencies in the limit 
of large noise. The stochastic coherence antiresonance, defined as a simultaneous occurrence of (i) 
the maximum of the coefficient of variation, and (ii) the minimum of the firing rate vs. the noise 
intensity, occurs over a wide range of parameter values, including monostable regions. Results of 
this work can be verified experimentally. 

PACS numbers: 87.19.1b,87.19.1c,87.19.1n 



I. INTRODUCTION 

The Hodgkin-Huxley (HH) model [Ij is a prime exam- 
ple of a resonant neuron. It was originally developed to 
explain experimental properties of the squid giant axon. 
Its behavior under the influence of constant, periodic 
Q and irregular external current (7| has been studied ex- 
tensively. It also served as a starting point for a number 
of reduced models [8l-[l0|. designed to preserve the key 
features, while being more amenable to large-scale nu- 
merical simulations. However, a full understanding of 
ah qualitative properties of its solutions has not been 
achieved yet. 

Deterministic HH equations can have both periodic 
and aperiodic, sometimes chaotic solutions [ll|. Theo- 
retical f]^^!^ and experimental f§, analysis revealed 
that near the excitation threshold two solutions, the fixed 
point and the limit cycle, may coexist. A simple such ex- 
ample is the HH model driven by a constant current. As 
the current magnitude is increased the neuron starts re- 
sponding at a preferred frequency, which is close to 50 
Hz for the original HH parameter set. 

The HH neuron in the presence of noise may display ei- 
ther resonant or antiresonant behavior, depending on the 
signal magnitude, frequency, and noise. The enhance- 
ment of weak signals by noise is known as stochastic res- 
onance p^|-[T7|. The opposite effect, stochastic antires- 
onance, in which a neuron's firing frequency is slowed 
down or even entirely stopped at some intermediate noise 
level, received some attention (l8l-,r20j recently. Depend- 
ing on model parameters this behavior is associated with 
bistability [18, "19] or multimodality of the response \2^\. 
A relation between bistability and antiresonance in inte- 
grator neurons was also studied recently (22| . 

In two previous papers on the HH model driven by pe- 
riodic sequence of short stimuli, we have shown that there 



is a transition between odd-only and all modes at high 
frequencies Ull [23| . This transition is located between 
the locked-in regions 2:1 and 3:1. The notation p : q 
means q output spikes for every p input current pulses. 
The edges of individual modes scale logarithmically in 
the vicinity of this singularity. This theoretical analysis 
agrees well with experimental data [24j. A natural ques- 
tion to ask is whether we can identify analogous parity 
transition, involving even-only modes on one side and all 
modes on the other side of the transition. A preliminary 
study [25] showed that the even and odd modes compete 
also near the main resonance. In the following we analyze 
the resonance regime in detail, looking for signatures of 
the even-all transition. 

Our aim is to map the response diagram of the HH 
neuron to a train of brief current pulses, rather than em- 
ulate typical in vivo scenarios. Stimuli in the form of brief 
pulses are better suited to reveal the internal dynamics 
of the neuron than signals varying on the time scale of 
the main resonance. In the case of the often used con- 
stant or sinusoidal input currents the neuron dynamics is 
obscured by the drive that is always different from zero. 
Gaussian noise is added to the model in order to reveal 
additional features of the model's dynamics. Although 
this particular type of randomness may not necessarily 
occur in a neuronal system we believe similar results will 
be obtained for any fast-varying irregular component of 
the stimulus. 



II. THE MODEL AND RESULTS 

We consider the stochastic HH model with the classic 
parameter set and rate constants 

= -iNa -Ik-Il+ lapp + C^{t), (1) 
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FIG. 1: Response diagram in the noiseless case for stimulation 
by rectangular current pulses of width r = 0.6ms and height 
Iq. The top (bottom) panel shows the high (intermediate) 
frequency regime. Bistable regions are marked "b". Some 
smaller bistable areas remain unlabeled. They are limited by 
a dotted line from above and a solid line from below. La- 
bels 1,2,3, and 4 mark regions with 1:1,2:1,3:1, and 4:1 mode 
locking, respectively. The resonance area near Ti = 17ms is 
dominated by the 3:1 state and higher order states. The sec- 
ond resonance near Ti = 34ms is dominated by the 2:1 mode. 
Higher modes appear near the threshold at Ti ~ 34ms. 



where noise is given by the Gaussian process < ^i{t) >= 
0, < ^it)^it') >= 2D6{t - t'), and D is expressed in 
mV^/ms. iNa, Ik, II, and lapp, are the sodium, potas- 
sium, leak, and external current, respectively. C = 
lliF/cm is the membrane capacitance. The input cur- 
rent is a periodic set of rectangular steps of height Iq 
and width 0.6ms, which is an order of magnitude be- 
low the resonant pulse width [20] and does not inter- 
fere with the neuron's internal dynamics. The deter- 
ministic (stochastic) HH equations are integrated us- 
ing the fourth-order Runge-Kutta algorithm (the second- 
order stochastic Runge-Kutta algorithm [26|), respec- 
tively. The simulations are carried out with the time step 
of 0.001ms and are run for 400s, discarding the initial 4s. 

Figure [1] shows the response diagram in the Ti — Iq 
plane without noise, where Ti is the stimulus period. The 
main resonance is located at T,- ~ 17ms. The second or- 



der resonance is present at Ti ~ 34ms and the third one 
at Ti ~ 51ms (not shown). The dotted line separates 
the region with a single solution from the area where 
two solutions coexist. Bistable solutions appear between 
the resonant regions near Ti = 11ms and Ti = 23ms. 
In these regimes the transition to excitability occurs via 
the subcritical Hopf bifurcation. The picture is more 
complicated at the frequencies above 250 Hz, where the 
quiescent state often coexists not with a limit cycle but 
with a set of irregular orbits. The boundaries of bistable 
regions are determined with a simple continuation algo- 
rithm. The initial conditions of each run with a new 
value of a bifurcation parameter are equal to the end 
values from the previous iteration. 

The firing rate fo/fi depends continuously on Iq be- 
tween the tip of the resonance at Ti ~ 17ms and the 
bistable area, which begins at Ti ~ 20ms. Here /o, fi 
is output, input frequency, respectively. For 14ms < 
Ti < 17ms bistable intervals with integer p/q ratio alter- 
nate with irregular response, whose long-time average is 
a continuous function of Iq and Ti. Earlier, Clay reported 
[13] an irregular graded response near the excitation edge 
within a revised HH model [28] stimulated by 1-ms rect- 
angular current pulses and trains of half-sine waves. This 
variability was attributed to the deterministic nonlinear 
dynamics of the model. 
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FIG. 2: The firing rate as a function of the current pulse 
height at the resonance Ti = 17.5ms. Here r = 0.6ms. Near 
the firing threshold fo is approximately a square root function 
of the pulse amplitude. Further away from the threshold the 
dependence of fo on Jo between the mode-locked states is 
irregular and nonmonotonic. 

For Ti > 1 7ms the firing rate is approximately a square 
root function of the deviation from the threshold current 
amplitude (see Fig. [2|). This dependence is characteris- 
tic of a saddle-node bifurcation [29]. The scaling fo/fi 
[Iq — IthY , with /3 = 1/2 is reminiscent of a mean-field 
second-order phase transition, with fd fi playing the role 
of an order parameter. The same exponent is obtained 
for a relaxation time near /e ~ 6.264/iA/cm^ for a con- 
stant current / where Ic is the value of / at the 
saddle-node bifurcation. Below In the neuron returns 
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FIG. 4: Sample voltage trace in the chaotic regime near the 
transition between even-only and all modes. The applied cur- 
rent pulses are shown with dotted lines. There is only one 
odd multiple (7:1) of the input period in this sample. Here 
T^ = 13.9ms and Iq = 10.75/xA/cm^. 



FIG. 3: Details of the response diagram at the resonance. 
Borders of even (odd) mode-locked states are shown with bro- 
ken (continuous) lines respectively. The dotted lines mark the 
transitions to repetitive firing. Areas enclosed by the dotted 
line and continuous or broken lines are bistable. Numbers 
2,. ..,5 indicate the states 2:1,...,5:1. The states beyond 7:1 
are not shown. Filled squares indicate the parity multimodal 
transition between even-only and all modes. 



to quiescence after emitting a series of spikes. The re- 
laxation time, defined as the time from the first to the 
last spike, diverges as (/c — /)^, where A ~ 1/2. Sim- 
ilar values of A were obtained for the Morris-Lecar and 
FitzHugh-Nagumo models [30]. 

Figure [3] shows the bifurcation diagram near the tip 
of the resonance. The synchronized states alternate with 
irregular firing. There are small bistable regions where 
the quiescent state coexists with a limit cycle. We have 
calculated their boundaries for states up to order 5:1. 
It is likely that they extend all the way to the tip of 
the resonance. The behavior in the intermediate regions 
is chaotic due to a competition between odd and even 
modes An example of V{t) dependence for a point 
located between states 2:1 and 3:1 is shown in Fig. [H 
This sample contains only one odd multiple of the input 
period. Odd interspike intervals vanish in the vicinity of 
the 2:1 state. Careful analysis of the interspike interval 
histograms (ISIH) reveals a transition between the set of 
even and the set of all modes. This parity transition of 
ISIHs in the HH model could be tested experimentally in 
a squid axon experiment, similar to the odd-all transition 
discovered earher |2i|, [23| • 

The transition between quiescence and chaotic firing 
marked by the dotted line in the intermediate zones be- 
tween the locked- in states occurs via period doubling (see 
Fig. [5|). The lowest and the highest values of V below 
Ti ~ 13.24ms belong to the 2:1 limit cycle. The line 
slightly below V = — 65mV, splitting above Ti = 13ms, 
is associated with the steady state. 

The even, odd multiples of Ti dominate near the states 
with even, odd p/q ratio, respectively (see Fig. [6l left). 
This is a general property of regimes of irregular firing. 
Starting from the 2:1 mode, the histogram weight is grad- 
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FIG. 5: The Poincare section of the membrane potential near 
the bistable region and irregular regime for current ampli- 
tude Iq = 10.65yU.A/cm^. The values of V{t) in each run are 
recorded at the points t = nTi, where n = 1, 2, .... The low- 
est and the highest flat branches on the left belong to the 
2:1 states. The middle one belongs to the steady state. It 
undergoes a period doubling bifurcation on approach to the 
irregular regime above Ti = 13.3ms. 

ually transferred to the 4:1 mode, which dominates for Ti 
slightly below the even- all transition point Tea- On the 
other side of Tea, = 1 over a narrow interval. 

The histogram of the dominant modes for Z) > is 
shown in Fig. [6] (right). Some weight is now transferred 
to lower modes. The 4:1 and 5:1 modes are still well 
pronounced over a range of Ti but they are now mixed 
with 2:1 and 3:1 modes, respectively. Noise extends the 
range of presence of high order modes around T^a, as 
expected. 

The average firing rate between nearby ml and (n+l):l 
states often has a narrow local minimum due to the pres- 
ence of slower modes. The minima are more pronounced 
near the excitation threshold along the left edge of the 
resonance tip in Fig. [3l where Ti < 17ms. One such 
valley is shown in Fig. [T] Here the 4:1 mode dominates 
in a narrow range of parameter values below the parity 
transition. 

The firing rate minimum is robust to small levels of 
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FIG. 6: Histograms of the lowest order modes for D = (left) and D = 10 ^ (right) at Iq = 10.75|uA/cm^. This is an irregular 
firing regime between the states 2:1 and 3:1. Data in the noisy case are Bezier-smoothed averages. 
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FIG. 7: Minimum of the firing rate between the 2:1 and 3:1 
states for the noiseless stimulus. The state 4:1 is located near 
the 2:1 mode. 
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FIG. 8: The firing rate as a function of Ti in the intermediate 
region between the 2:1 and 3:1 states for D = 10~^ and Ti = 
10.75ms. The competition of odd and even modes manifests 
itself as a minimum of the firing rate. 



noise, as can be seen in Fig. [HI The minimum of fo 
occurs for Ti > T^a, where modes of both parity are 
available. In the odd-all transition of the high-frequency 
regime [23| the minimum fo was similarly located close 
to the transition, on the side of all modes. 

It is well known that the HH neuron has a tendency to 
spike in bursts when subjected to a noisy stimulus in a 



bistable regime. If the deterministic system is prepared in 
one of high order bistable states from Fig. [Sj the addition 
of noise results in slow bursts, where the ISI within a 
burst is a high integer multiple of T^. An example is 
shown in Fig. ^ where each burst's ISI is equal to ST^. 
We expect slow bursts of order 6:1 and higher to be found 
in a more detailed calculation. 

For stimuli below the deterministic threshold noise ac- 
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FIG. 9: Sample voltage trace in the bistable regime of the 
state 3:1 under the influence of small noise. The most common 
form of response consists of 3:1 bursts separated by longer 
silent intervals. Here Ti = 13.9ms, Iq = 10.75/xA/cm^, and 
D = 10"^ 
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FIG. 10: The firing rate as a function of Ti for different noise 
levels and / = 10/iA/cm^, about 0.1/xA/cm^ below the deter- 
ministic threshold. 



tivates firing with maximum response at the resonant fre- 
quency (see Fig. [TU]) . When noise is larger, excitations 
into the 2:1 state occur more frequently and the maxi- 
mum of the firing rate shifts toward Ti = 20ms. This 
maximum is a result of the interplay between the 1:1 
mode and the 2:1 mode. The relative frequency of par- 
ticipation of the 2:1 mode grows sharply for intermediate 
noise levels in the neighborhood of the bistable regime 
for T, ~ 20ms. 

For signals with a periodic subthreshold component 
noise may play the role of a frequency selector. In a net- 
work with some elements firing in unison at the resonant 
frequency the remaining uncorrelated neurons also play 
an important role. The intensity of their background ac- 
tivity may select the firing rate of a network. 

The dependence of Cy on Ti for fixed noise intensity 
is shown in Fig. [TTl As expected, the most coherent 
response occurs near the resonance. The minimum near 
Ti = 20ms for suprathreshold inputs follows the left edge 
of the 1:1 region from Fig. [T] For stimuli below threshold 
there is a deep minimum at the resonance and a large 
maximum in the bistable regime for Ti > 20ms. 

For suprathreshold stimulation in the resonance regime 
the firing rate in general increases monotonically as a 
function of D with the exception of states having p/q 
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FIG. 11: The coefficient of variation as a function of Ti for 
different pulse amplitudes. The noise intensity if fixed at D = 
0.2. The local minimum at larger Jo follows the left edge of 
the 1:1 regime. 



values close to 1 (see Fig. [TT]). Starting in the 1:1 state 
at D = 0, increasing noise slows down the response by 
annihilating some of the action potential spikes. The 
minimum fo/fi is reached at an intermediate D. In the 
limit of large noise fo approaches the inverse of the refrac- 
tory period. The Cy maxima are caused by redistribu- 
tion of histogram weight among several principal modes. 
This effect was described earlier as stochastic coherence 
antiresonance [2l|. Maximization of spike train incoher- 
ence was shown earlier to occur in the FitzHugh-Nagumo 
[31] and leaky integrate-and-fire [32] models. However in 
those studies the maxima of Cy were not accompanied 
by the minima of fo/fi- 

Figure [13] shows the response diagram in the pres- 
ence of noise. The tip of the resonance is broadened 
and shifted to higher frequencies. The 3:1 state vanishes 
completely for D slightly larger than 10~^. Higher-order 
states are more sensitive to noise and are quickly washed 
out. The bistable area at the edge of the 2:1 state shrinks 
more gradually. When bistability is finally eliminated, 
fo remains discontinuous over part of that border. For 
D = 10~^ the discontinuity occurs below Ti ~ 10.5ms. 
Above Ti ~ 10.5ms the firing rate is a continuous func- 
tion of Iq. The behavior in the immediate vicinity of 
Ti = 10.5ms seems to be weakly irregular. A more de- 
tailed study would be needed for a proper description of 
this area. The loss of stability by the 2:1 state in favor 
of the quiescent state and further transition to bursting 
are illustrated in Fig. [TH 

The reaction to noise at high frequencies is qualita- 
tively different from behavior at near-resonant frequen- 
cies (see Fig. [T5|) . In the top panel of Fig. [15] the firing 
rate drops almost everywhere for small D. This is easy to 
understand if we recall that the odd-all multimodal tran- 
sition [U, [23I occurs just below Iq = 20/iA/cm^. The 
slow modes are easily available in this regime and small 
perturbations suffice to switch the system among differ- 
ent trajectories. The firing rate drops over the entire 3:1 
plateau due to the appearance of the 5:1, 7:1 and other 
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FIG. 12: The firing rate and Cy as a function of D for several 
values of /o and Ti = 17ms. Maximum of Cv at intermediate 
values of D is due to the multimodal distribution of interspike 
intervals. 



odd modes. For larger D also the even modes are sam- 
pled by the system and fo/fi increases. 

At moderate, near-resonance frequencies, the central 
part of the 2:1 plateau in the lower part of Fig. [15] is 
more robust to noise. The average fo is preserved for D 
up to 0.1. Noise induces both the 1:1 mode as well as the 
3:1 and slower modes. However the participation rate of 
the 1:1 mode is balanced by the slower modes in such a 
way that average fo/fi remains close to 2. This type of 
resilience to noise, characteristic of the resonance regime, 
can also be seen in other mode-locked states. 

Finally, let us briefly describe the role played by the 
current pulse width r. Figure [16] shows the excita- 
tion threshold as a function of r for a resonant drive, 
Ti = 17ms. Initially the threshold decreases as the in- 
verse of T. For T > 8ms the tip of the resonance (shown 
in Fig. [3]) gradually becomes narrower and vanishes. 
The threshold rises again and reaches maximum near 
T = 17ms, when I{t) = const. Bistability appears above 
T ~ 16ms. It is clear that stimulation by a constant cur- 
rent forces the neuron into an antiresonant regime. The 
bifurcation diagram in Fig. [16] implies that the addi- 
tion of any charge-unbalanced component to the constant 
stimulus takes the system away from this antiresonant 
limit. The precise nature of such a component, whether 
deterministic or stochastic, is not important. What mat- 
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FIG. 13: The response diagram in the presence of noise. Solid 
lines are the boundaries of the main mode-locked states with- 
out noise (see also Fig. [T]). The broken and dotted lines are 
boundaries for D = 10~^ and 10~^, respectively. 
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FIG. 14: Sample dependence of V{t) for three intensities of 
noise, from top to bottom: D = 0.001,0.01, and 0.03. The 
spiking action is switched off for intermediate values of D. 
Here Ti = 11ms, Iq = lO.S/xA/cm^- 



ters is the amount of charge delivered within approxi- 
mately 4ms from the stimulus onset. 



III. CONCLUSIONS 

We studied the response of the HH neuron to a periodic 
pulse current with a Gaussian noise. The global bifurca- 
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FIG. 15: Typical behavior of the average firing rate in the 
regime of high (top) and moderate frequency (bottom) for 
different noise levels. At high frequencies and small noise 
the response is initially slowed down over an entire plateau, 
starting from both edges. At the lower diagram /o is held 
steady in the central part of the plateau, decreasing at the 
left edge and increasing at the right edge. The average firing 
rate in the central part of the plateau is preserved up to D 2± 
0.1. The irregular nonmonotonic behavior in the deterministic 
case is smoothed out for very small levels of noise. 

tion diagram in the Ti — Iq plane has a rich structure 
near the excitation threshold, where resonant regimes al- 
ternate with antiresonant ones. The model is bistable be- 
tween the resonances and in the limit of high-frequency 
stimulation. 

The firing rate is a continuous function of Iq for 
Tres ^ Ti < 20ms. The scaling of the firing rate with 
{Iq — IthY^'^ is a signature of a saddle- node bifurcation 
at the threshold. For Ti < Tres bistable regions are sepa- 
rated by areas of irregular response with no well defined 
threshold and approximately continuous dependence of 
fo/ fi on input parameters. As Ti approaches Tres from 
below, the subcritical Hopf bifurcation gradually softens. 
Bistable regions occupy smaller portions of the parame- 
ter space and disappear before Ti — Tres- The change 
of the type of neuronal excitability is important in the 
context of preventing the so called dynamical diseases, 
such as epilepsy, Alzheimer's, and Parkinson's disease. 
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FIG. 16: Bifurcation diagram near the excitation threshold as 
a function of the current pulse width. The input period is set 
at Ti = 17ms. At r = 17ms the input current is constant. The 
details of the antiresonant limit with bistable behavior are 
shown in the inset. Solid lines mark the excitation threshold. 
The dotted lines are the boundaries of the 1:1 state. 

It was shown in the HH [sH, [s^] and the Hindmarsh-Rose 
(SS*] models that this could be obtained by introducing an 
additional control function. We demonstrated that the 
threshold behavior can also be altered by selecting the 
frequency and amplitude of external current. 

For subthreshold stimuli noise enables spiking in the 
vicinity of the resonance. The stimulus frequency, for 
which the maximum firing rate is obtained, depends on 
the magnitude of noise. For large D, the maximum fo oc- 
curs for Ti > Tres- It would be interesting to investigate 
the same regime in a HH network where the connectivity 
pattern of neurons as well as their individual properties 
might lead to the emergence of subpopulations of neurons 
firing with different average frequencies in response to a 
correlated input with background noise. It was shown 
earlier by several authors that both the background noisy 
activity [36, 37j and correlated inputs [38j are important 
in explaining the neuronal response in vivo. Qualitative 
results of this work do not depend on the precise func- 
tional form of the current pulse provided the width of 
each pulse does not exceed 4 ms. This invariance of the 
bifurcation diagram in the Ti — Iq plane is not surpris- 
ing because for short pulses the HH neuron's threshold is 
determined by the amount of charge delivered per pulse 
[39]. 

We have also found a new even-all multimodal transi- 
tion occurring between the states 2:1 and 3:1 close to the 
main resonance. For input period Ti below this singular- 
ity only even response modes exist. Both even and odd 
modes appear above the transition. This effect is accom- 
panied by a minimum of the firing rate, located close to 
the parity transition, on the side of all modes. Similar 
transitions may exist in other excitable systems. 

Our results are also relevant to studies of auditory 
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nerve fiber responses to electric stimulation (40l-l43|. At 
low stimulation rates auditory nerve fibers fire regularly 
and are locked in to applied stimulus. At high stimulation 
rates these fibers respond irregularly. Some researchers 
attributed this effect to physiological noise [44, 45] . How- 
ever, an analysis within the FitzHugh-Nagumo model 
showed that the firing irregularities at high frequencies 
may be caused by deterministic dynamical instability 
[421 . |43|. We have found similar instabilities in the HH 
model. They appear at high stimulation frequencies and 



along the excitation threshold. Both of these regimes are 
relevant to studies of hearing sensitivity. 
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